﻿
class caja
{
private:

public:
	double lbox[3];//Lado de la caja
	particula *part; //Array de particulas
	int npart;//Numero de particulas
	double phi;//Concentracion en volumen
	double rho;//Concentracion en masa
	double temp;//Temperatura de caja
	double vol;//Volumen
	double rcut;//Distancia de corte de interaccion
	///Para saber tanto por ciento de movimientos aceptados
	double ntry;//Numero de intentos
	double nok; //Movimientos aceptados

	///Parametros lara la link-cell (vecinos)
	int ncell;//Numero de celdas 
	int ncell2;//ncell^2
	int ncell3;//ncell^3
	double lgrid;//Lado de celda
	/*head[i] contiene el primer indice de particula contenido en celda i
	head[i]=-1 si no hay particulas en celda i*/
	int *head;
	//Contiene la lista de particulas enlazadas
	int *link;
	/*Ejemplo:
	celda i
	p1=head[i]
		- p1 es el indice de la primera particula en celda i. 
		- si no hay particulas en celda i, p1=-1
	p2 = link[p1]
		- p2 es el indice de la segunda particula en celda i
		- p2=-1 si no hay mas particulas en la celda
	etc...
	Por lo tanto, en la celda i habrá
	[p1,p2,..., pk] particulas	
	*/
	//Energias
	double epoten;//Energia potencia por particula en la caja
	double virial;
	double ekin;//Energia cinetica por particula en la caja
	//Control de la velocidad (termostato)
	int nacc;///Paso de la simulacion
	double delta;//Paso de longitud

	caja()
	{

	}

	caja(int np, double xT, double xphi, double xrcut, double del)
	{
		///Empleamos paso de tiempo Dt para generar las posiciones posOld
		///de las particulas empleando las velocidades
		ntry=0;
		nok=0;
		delta=del;
		npart=np;
		temp =xT;
		phi=xphi;
		rcut = xrcut;
		rho = phi*6.0/pi;
		double bx = pow(3.1416*npart/(6.*phi),1./3.);

		for(int i=0;i<3;i++)
			lbox[i]= bx;
		part = new particula[npart];
		link = new int[npart];
		SetMeshPart();
		vol=lbox[0]*lbox[1]*lbox[2];

		ncell=int(lbox[0]/rcut);
		ncell2=ncell*ncell;
		ncell3=ncell*ncell*ncell;
		lgrid= lbox[0]/double(ncell);
		cout<<"Numero de celdas: "<<ncell<<endl;
		cout<<"Lado de celda: "<<lgrid<<endl;
		head = new int[ncell*ncell*ncell];
	}

	///Reinicia la simulacion desde un archivo
	caja(char filename[50], double del)
	{
		ntry=0;
		nok=0;
		delta=del;
		double x;
		///Empleamos paso de tiempo Dt para generar las posiciones posOld
		///de las particulas empleando las velocidades
		ifstream MiArchivo (filename,ios::in);
		if (MiArchivo.is_open())
		{
			///Cargamos las caracteristicas de la caja
			MiArchivo>>npart>>phi>>temp>>rcut;
			rho = phi*6.0/3.1416;
			double bx = pow(3.1416*npart/(6.*phi),1./3.);
			for(int i=0;i<3;i++)
				lbox[i]= bx;
			part = new particula[npart];
			link = new int[npart];
			SetMeshPart();
			vol=lbox[0]*lbox[1]*lbox[2];

			ncell=int(lbox[0]/rcut);
			ncell2=ncell*ncell;
			ncell3=ncell*ncell*ncell;
			lgrid= lbox[0]/double(ncell);
			cout<<"Numero de celdas: "<<ncell<<endl;
			cout<<"Lado de celda: "<<lgrid<<endl;
			head = new int[ncell*ncell*ncell];
			///Cargamos las posiciones de las particulas
			for(int i=0;i<npart;i++)
			{
				for(int j=0;j<3;j++)
				{
					MiArchivo>>part[i].pos[j];
					//Posiciones anteriores de las particulas
					part[i].posOld[j]=part[i].pos[j];
				}
			}
			MiArchivo.close();
		}

	}

	~caja()
	{
		delete [] part;
	}

	///Imprime la configuracion
	void printConfig(void);
	///Coloca las particulas de forma aleatoria en una caja
	void SetRandomPart(void);
	///Coloca las particulas en red
	void SetMeshPart(void);
	///Imprime la informacion de la particula considerada
	void followPart(int ind);
	///Calcula la linkCell
	void linkCell(void);
	///Calcula la lista de vecinos
	void listaVecinos(void);
	///Aplica condiciones contorno periodicas a particula ind
	void boundary(int ind);
	///Calcula la distancia cuadratica entre particulas
	double distancia2(int ind1,int ind2);
	///Calcula la distancia entre particulas
	double distancia(int ind1,int ind2);
	///Imprime el archivo XYZ de la configuracion de particulas
	void printXYZfile(char namefile[30], int i);
	//Selecciona una particula al azar y la desplaza
	void randomMovement(void);
	//Inicializa la energia de cada particula
	void calcEnergia(void);
	//Calcula la energia de la particula ind
	double calcEnergiaParticle(int ind);
	///Elimina a una particula de una celda
	void removeParticleFromCell(int ind);
	///Inserta particula en una celda
	void insertParticleInCell(int icell, int ind);
	int testLinkCell(void);

};

int caja::testLinkCell(void)
{
	int cual, nnn;
	nnn=0;
	for(int i=0;i<ncell3;i++)
	{
		cual=head[i];
		while(cual!=-1)
		{
			nnn++;
			cual=link[cual];
		}
	}
	return nnn;
}

void caja::randomMovement(void)
{
	double alea = (double(rand())/RAND_MAX);
	int index = int(alea*npart);
	if (index==npart) index--;
	double posR;
	int mc[3];
	int icellOld;
	int icellNew;
	double energia;

	//Actualizamos el numero de intentos
	ntry++;
	//Celda en la que se encontraba la particula
	icellOld=part[index].celda;
	//Desplazamos aleatoriamente la particula
	for (int j=0;j<3;j++)
	{
		posR=(part[index].pos[j])+delta*(double(rand())/RAND_MAX-0.5);
		part[index].posOld[j]=part[index].pos[j];
		part[index].pos[j]=posR;
	}
	//Condiciones de contorno periodicas
	boundary(index);
	//Localizamos la celda en la que se encuentra la particula
	for(int j=0;j<3;j++)
	{
		mc[j] = int((part[index].pos[j]+0.5*lbox[0])/lgrid);
	}
	icellNew = int(mc[0]*ncell2+mc[1]*ncell+mc[2]);
	
	if(icellNew!=icellOld)
	{
		//Actualizamos las coordenadas de la celda donde se encuentra la particula
		for(int j=0;j<3;j++)
		{
			part[index].celdaxyz[j]=mc[j];
		}
		//La particula ha cambiado de celda
		//Sacamos a la particula de la celda en la que se encontraba
		removeParticleFromCell(index);
		//Introducimos a la particula en la nueva celda
		insertParticleInCell(icellNew, index);
		//Acualizamos el objeto 'particula' con la nueva celda
		part[index].celda=icellNew;		
	}

	///Calculamos la nueva energia de la particula
	energia = calcEnergiaParticle(index);
	if((double(rand())/RAND_MAX)<exp(-(energia-part[index].energia)/temp))
	{
		//Actualizamos el numero de intentos satisfactorios
		nok++;
		//cout<<"aceptado"<<endl;
		//Aceptamos el movimiento y actualizamos la energia
		part[index].energiaOld=part[index].energia;
		part[index].energia=energia;
	}
	else
	{
		//cout<<"rechazado"<<endl;
		//Rechazamos el movimiento y colocamos a la particula en su posicion anterior
		for(int j=0;j<3;j++)
		{
			part[index].pos[j]=part[index].posOld[j];
		}
	}
}

void caja::printConfig(void)
{
	string posPart="";
	char aux[100];
	for(int i=0;i<npart;i++)
	{
		sprintf(aux,"%f%s%f%s%f%s",part[i].pos[0]," ",part[i].pos[1]," ",part[i].pos[2],"\n");
		posPart+=aux;
	}
	sprintf(aux,"%s%d%s","conf_",nacc,".txt");
	ofstream MiArchivo (aux,ios::out);
	if (MiArchivo.is_open())
	{
		MiArchivo<<npart<<" "<<phi<<" "<<temp<<" "<<rcut<<endl;
		MiArchivo<<posPart;
		MiArchivo.close();
	}
}

double caja::calcEnergiaParticle(int index)
{
	int c, c1;
	int iii,jjj;
	int shift[3];
	int ank, ank2, ankh;
	ank=6;
	ank2=2*ank;
	ankh=ank/2;
	double dist2,rqInv,tsn, potdm, ePoten;
	double conrep= 4.0*pow(1.0,ank2); //Parte repulsiva del potencial Leonard-Jones
	double conat= 4.0*pow(1.0,ank); //Parte atractiva del potencial Leonard-Jones
	double rcut2=rcut*rcut;

	ePoten=0.;

	c = part[index].celda;
	///Recorremos las particulas de esta celda y las celdas vecinas
	for (int ixred=part[index].celdaxyz[0]-1;ixred<=part[index].celdaxyz[0]+1;ixred++)
	{
		/*shift: condiciones de contorno periodicas tenidas en cuenta para
		localizar las celdas vecinas*/
		if (ixred < 0)
		{
			shift[0] = int(ixred+ncell);
		}
		else
		{
			if (ixred>=ncell)
			{
				shift[0] = int(ixred-ncell); 
			}
			else
			{
				shift[0] = int(ixred);
			}
		}
		for (int iyred=part[index].celdaxyz[1]-1;iyred<=part[index].celdaxyz[1]+1;iyred++)
		{
			if (iyred < 0)
			{
				shift[1] = int(iyred+ncell);
			}
			else
			{
				if(iyred>=ncell)
				{
					shift[1] = int(iyred-ncell);
				}
				else
				{
					shift[1] = int(iyred);
				}
			}
			for (int izred=part[index].celdaxyz[2]-1;izred<=part[index].celdaxyz[2]+1;izred++)
			{
				if (izred < 0)
				{
					shift[2] = int(izred+ncell);
				}
				else
				{
					if(izred>=ncell)
					{
						shift[2] = int(izred-ncell);
					}
					else
					{
						shift[2] = int(izred) ;
					}
				}
				///Calculamos el indice c1 de la celda vecina
				c1 = int(((shift[0]+ncell)%ncell)*ncell2+((shift[1]+ncell)%ncell)*ncell+((shift[2]+ncell)%ncell));
				///El indice index ahora es iii
				iii = index;//bx->head[c];

				///Buscamos la primera particula jjj de la celda c1
				jjj =head[c1];


				while (jjj != -1) ///jjj=-1->No hay particulas en la celda
				{
					if(iii!=jjj)
					{

						///Calculamos la distancia con condiciones periodicas entre iii y jjj
						dist2= distancia2(iii,jjj);
						//Si la distancia es menor que el cutoff calculamos la interaccion
						if (dist2 < rcut2 )
						{
							///Calculamos la fuerza entre iii y jjj
							rqInv =1.0/dist2;//1/dist^2
							tsn=pow(rqInv,ankh);
							potdm = tsn*(conrep*tsn-conat);
							//Acualizamos la energia potencial
							ePoten+=potdm;
						}
					}
						//Pasamos a la siguiente particula de la celda c1
						jjj = link[jjj];
				}
			}
		}
	}
	return ePoten;
}

void caja::calcEnergia(void)
{
	int c, c1;
	int iii,jjj;
	int shift[3];
	int ank, ank2, ankh;
	ank=6;
	ank2=2*ank;
	ankh=ank/2;
	double dist2,rqInv, tsn, potdm, ePoten;
	double conrep= 4.0*pow(1.0,ank2); //Parte repulsiva del potencial Leonard-Jones
	double conat= 4.0*pow(1.0,ank); //Parte atractiva del potencial Leonard-Jones
	double rcut2=rcut*rcut;
	//Correccion al cutoff del potencial
	double LJcut= pow(1.0/(rcut*rcut),ankh)*(pow(1.0/(rcut*rcut),ankh)*conrep-conat);

	for (int i=0;i<npart;i++)
	{
		part[i].energia=0;
		part[i].energiaOld=0;
	}

	ePoten=0.;
	//Recorremos todas las celdas
	for (int ix=0;ix<ncell;ix++)
	{
		for (int iy=0;iy<ncell;iy++)
		{
			for (int iz=0;iz<ncell;iz++)
			{
				///Calculamos el indice de la celda (ix,iy,iz)
				c = int(ix*ncell2+iy*ncell+iz);
				///Recorremos las particulas de esta celda y las celdas vecinas
				for (int ixred=ix-1;ixred<=ix+1;ixred++)
				{
					/*shift: condiciones de contorno periodicas tenidas en cuenta para
					localizar las celdas vecinas*/
					if (ixred < 0)
					{
						shift[0] = int(ixred+ncell);
					}
					else
					{
						if (ixred>=ncell)
						{
							shift[0] = int(ixred-ncell); 
						}
						else
						{
							shift[0] = int(ixred);
						}
					}
					for (int iyred=iy-1;iyred<=iy+1;iyred++)
					{
						if (iyred < 0)
						{
							shift[1] = int(iyred+ncell);
						}
						else
						{
							if(iyred>=ncell)
							{
								shift[1] = int(iyred-ncell);
							}
							else
							{
								shift[1] = int(iyred);
							}
						}
						for (int izred=iz-1;izred<=iz+1;izred++)
						{
							if (izred < 0)
							{
								shift[2] = int(izred+ncell);
							}
							else
							{
								if(izred>=ncell)
								{
									shift[2] = int(izred-ncell);
								}
								else
								{
									shift[2] = int(izred) ;
								}
							}
							///Calculamos el indice c1 de la celda vecina
							c1 = int(((shift[0]+ncell)%ncell)*ncell2+((shift[1]+ncell)%ncell)*ncell+((shift[2]+ncell)%ncell));
							///Buscamos la primera particula iii de la celda c
							iii = head[c];
							while (iii != -1) ///iii=-1 ->No hay particulas en la celda
							{
								///Buscamos la primera particula jjj de la celda c1
								jjj =head[c1];
								while (jjj != -1) ///jjj=-1->No hay particulas en la celda
								{
									//Evitamos el conteo doble de las interacciones
									if (iii < jjj) 
									{
										///Calculamos la distancia con condiciones periodicas entre iii y jjj
										dist2= distancia2(iii,jjj);
										//Si la distancia es menor que el cutoff calculamos la interaccion
										if (dist2 < rcut2)
										{
											///Calculamos la fuerza entre iii y jjj
											rqInv =1.0/dist2;//1/dist^2
											tsn=pow(rqInv,ankh);
											potdm = tsn*(conrep*tsn-conat);
											//Acualizamos la energia potencial
											ePoten+=potdm;
							/*				if(potdm>100)
											{
												cout<<"ss"<<endl;
											}*/
											//Actualizamos la energia individual de las particulas
											if(iii==5 || jjj==5)
											{
												//cout<<"aaa"<<endl;
											}
											part[iii].energia+=potdm;
											part[jjj].energia+=potdm;
											
										}
									}
									//Pasamos a la siguiente particula de la celda c1
									jjj = link[jjj];
								}
								//Pasamos a la siguiente particula de la celda c
								iii = link[iii];
							}
						}
					}
				}
			}
		}
	}
	//Introducimos la correccion del cutoff del potencial
	ePoten-=LJcut*(npart)*((npart)-1)/2.0;
	epoten = ePoten/npart;
}

void caja::SetRandomPart(void)
{
	double pos[3];

	for(int i=0;i< npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			pos[j]=lbox[j]*double(rand())/RAND_MAX;
		}
		part[i].setPosition(pos);
		boundary(i);
	}
}

void caja::SetMeshPart(void)
{
	int np = int(ceil(pow(npart,1./3.)));
	double lx = lbox[0]/double(np);
	int n=0;
	double pos[3];
	double drift=0.5*lx;

	pos[0]=-0.5*lbox[0]+drift;
	for(int ix=0;ix<np;ix++)
	{
		pos[1]=-0.50*lbox[0]+drift;
		for(int iy=0;iy<np;iy++)
		{
			pos[2]=-0.50*lbox[0]+drift;
			for(int iz=0;iz<np;iz++)
			{
				for(int k=0;k<3;k++)
				{
					pos[k]+=(double(rand())/RAND_MAX-0.5)*0.01;
				}
				part[n].setPosition(pos);

				boundary(n);
				pos[2]+=lx;
				n++;

			}
			pos[1]+=lx;
		}
		pos[0]+=lx;
	}
	cout<<"lbox: "<<lbox[0]<<endl;
	cout<<"lx: "<<lx<<endl;
	cout<<"np: "<<np<<"  "<<int(ceil(pow(npart,1./3.)))<<endl;
}

void caja::followPart(int ind)
{
	cout<<"pos: ";
	for(int i=0;i<3;i++)
		cout<<part[ind].pos[i]<<"  ";
	cout<<endl;
}

void caja::linkCell(void)
{
	int mc[3], c;
	//Reset the headers, head
	for (int i=0;i<ncell3;i++)
	{
		head[i] = -1;
	}

	///Scan atoms to construct headers, head, & linked lists, lscl
	for (int i=0;i<npart;i++)
	{
		///Vector cell index to which this atom belongs
		for (int j=0;j<3;j++)
		{
			mc[j] = int((part[i].pos[j]+0.5*lbox[0])/lgrid);
			//Coordenadas de la celda en la que se encuentra la particula
			part[i].celdaxyz[j]=mc[j];
		}
		///Translate the vector cell index, mc, to a scalar cell index"""
		c = int(mc[0]*ncell2+mc[1]*ncell+mc[2]);
		
		
		//Introducimos la celda en la que se encuentra la particula en el objeto 'particula'
		part[i].celda=c;
		///Link to the previous occupant (or EMPTY if you're the 1st)
		link[i] = head[c];
		///The last one goes to the header"""
		head[c] = i;
	}

	int ntotal=0;
	for(int i=0;i<ncell3;i++)
	{
		int cual = head[i];
		int n=0;
		while(cual!=-1)
		{
			cual=link[cual];
			n++;
		}
		ntotal+=n;
	}
	cout<<"Numero total de particulas colocadas en Mesh "<<ntotal<<endl;
}

void caja::removeParticleFromCell(int ind)
{
	int icell = part[ind].celda;
	int i = head[icell];
	int j, iaux;
	if(i==ind)
	{
		j=link[i];
		head[icell]=j;
	}
	else
	{
		while (i!=-1)
		{
			iaux = link[i];
			if (iaux==ind)
			{j=link[iaux];}
			else
			{j=iaux;}
			link[i]=j;
			i=j;
		}
	}

	int nnn = testLinkCell();
	if(nnn!=npart-1)
	{
		cout<<"aaa"<<endl;
	}

}

void caja::insertParticleInCell(int icell, int ind)
{
	int nnn = testLinkCell();
	int iold;
	int i = head[icell];
	if (i==-1)
	{
		head[icell]=ind;
		link[ind]=-1;
	}
	else
	{
		while (i!=-1)
		{
			iold=i;
			i = link[i];
			if (i==-1)
			{
				link[iold]=ind;
				link[ind]=-1;
				i=-1;
			}
		}
	}

	nnn = testLinkCell();
	if(nnn!=npart)
	{
		cout<<"aaa"<<endl;
	}
}

void caja::listaVecinos(void)
{
	double r2, rc2;
	double shield = 0.5;
	rc2=(rcut+shield)*(rcut+shield);
	for (int i=0;i<npart;i++)
	{
		part[i].nvecinos=0;
	}

	for (int i=0;i<npart-1;i++)
	{
		for(int j=i+1;j<npart;j++)
		{
			r2 = distancia(i,j);
			if(r2<rc2)
			{
				part[i].lstVecinos[part[i].nvecinos]=j;
				part[i].nvecinos++;
				part[j].lstVecinos[part[j].nvecinos]=i;
				part[j].nvecinos++;


			}
		}
	}
}

void caja::boundary(int ind)
{
	for(int i=0;i<3;i++)
	{
		int nper = int(dabs(part[ind].pos[i])/(0.5*lbox[i])); //si esta dentro de la caja 0; si esta fuera 1
		double anper = sign((double(nper)*lbox[i]),part[ind].pos[i]); //0 dentro de caja, 1 � -1 si esta fuera
		//Condiciones de contorno peri�dicas
		part[ind].pos[i]-=anper;
	}
}

double caja::distancia2(int ind1,int ind2)
{
	//Calcula la distancia entre las particulas i y j
	double dr[3];
	double dis=0;
	for (int i=0; i<3;i++)
	{
		dr[i] = part[ind1].pos[i]-part[ind2].pos[i];
		///Condiciones de contorno periodicas
		dr[i]-=int(dr[i]/(0.5*lbox[i]))*lbox[i];
		dis+=dr[i]*dr[i];
		//cout<<dr[i]<<" ";
	}
	//cout<<dis<<endl;
	return	 dis;
}

double caja::distancia(int ind1,int ind2)
{
	return	 sqrt(distancia2( ind1,ind2));
}

void caja::printXYZfile(char namefile[30], int i)
{
	char name[50];
	sprintf (name, "%s%s", namefile, ".xyz");
	ofstream MiArchivo (name,ios::out|ios::app);
	if (MiArchivo.is_open())
	{
		MiArchivo << npart<<endl;
		MiArchivo << lbox[0]<<"  "<<i<<endl;
		for(int i=0;i<npart;i++)
		{
			MiArchivo<<"O "<<part[i].pos[0]<<"  "<<part[i].pos[1]<<"  "<<part[i].pos[2]<<endl;
		}
		MiArchivo.close();
	}
}

